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PREFACE 


This  research  of  the  three-dimens  ipfial  viscous  vortex/wing 
interaction  theory  was  performed  under  Contract  No.  NOOOI^-7^- 
C-0151  for  the  Office  of  Naval  Research.  The  work  was  done  in 
the  Lockheed-Georg ia  Company  Advanced  Flight  Sciences 
Department,  managed  by  Mr.  H.  R.  Leslie.  Major  contributors  to 
the  effort  came  from  consultants  Drs.  J.  F.  Nash  and 
R.  M.  Scruggs  of  Sybucon,  Inc.  The  program  monitor  is 
Dr.  R.  E.  Whithead,  Office  of  Naval  Research,  Arlington,  Va. 
This  work  is  continuing  currently. 


SUMMARY 


A computational  fluid  mechanical  model  developed  for  analyzing  the  loads 
and  flow  mechanism  due  to  vortex/wing  interactions  has  been  further  improved 
and  evaluated.  The  model  provides  vorticity,  velocity  vectors,  and  pressures 
for  viscous'  flow  on  the  upper  surface  of  a wing  with  a leading  edge  vortex. 
Viscous  parametric  investigations  have  been  conducted  showing  the  effects  of 
such  parameters  as  leading  edge  sweep,  leading  edge  feeding  vorticity,  and 
viscosity.  Effects  of  grid  size  and  various  boundary  conditions  are 
presented. 

Boundary  conditions  are  determined  from  potential  flow  models  by  an 
iteration  procedure.  Directions  and  an  example  of  the  interaction  procedure 
is  presented.  A simple  65°  delta  wing  at  22°  angle  of  attack  is  used  for  the 
analyses.  Iterations  for  this  model  are  not  complete,  but  the  results  are 
instructive  for  analyzing  other  models  and  completing  the  65°  wing  investi- 
gation. The  method  is  not  limited  to  thin  delta  wings.  Arbitrary  planforms 
can  be  evaluated.  Thick  wings  with  rounded  leading  edges  are  under  current 
investigation. 


Experimental  laser  velocimeter  data  are  presented  for  boundary 
conditions  at  one  station  on  the  65°  delta  wing  at  22°  angle  of  attack. 
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1.  INTRODUCTION 

The  favorable  lifting  effects  of  leading  edge  vortices  on  delta  wings 
have  been  under  investigation  for  many  years.  More  recently,  other  types 

'j 

of  favorably  interfering  vortices  have  been  noteti.  These  consist  of  stream- 
wise  vortices  formed  from  leading  edge  "strakes",  outboard  leading  edge 
extensions  ("snags")  as  on  the  F-4  (aircraft)  and  trailing  vortices  from 
canard  configurations.  Now  that  interfering  vortices  are  recognized  aids, 
as  well  as  detriments,  to  lift  and  stability  augmentation,  considerable 
emphasis  is  being  placed  on  gaining  knowledge  concerning  their  formation 
and  stabi 1 ity. 

There  have  been  many  attempts  to  model  the  flow  for  delta  wings  with 
leading  edge  vortices.  Matoi^  has  provided  a good  survey  of  many  of  these 
attempts.  The  efforts  of  most  have  been  limited  to  potential  flow  models, 
assuming  rotational  elements  for  viscous  effects  such  as  the  leading  edge 
feeding  sheet.  The  basic  assumption  of  slender-body,  conical  flow  have  been 
followed  by  most  investigators  such  as  Mangier  and  Smith^.  Results  of  these 
early  investigator  e laid  a good  foundation  for  further  research,  but 
they  are  not  su  in  obtaining  forces  and  moments  over  a general  range 

of  configurat  >namus3,  with  the  leading  edge  suction  analogy,  has 

good  success  in  predicting  vortex  lift.  His  method  does  not  describe  the 
pressures  or  flow  mechanism.  Good  results  were  obtained  most  recently  by 
Weber,  et  al.**  and  Rao  and  Nathnan^,  where  the  trailing  edge  conditions  of 
delta  wings  are  included  in  the  boundary  conditions.  Kandil^  has  also 
provided  a good  potential  flow  model  and  convergence  criteria.  The  latter 
have  limits  that  restrain  their  application  to  certain  configurations  as  well 
as  raise  doubt  concerning  their  accuracy.  These  limits  are  given  below. 

(1)  Thin  wings  only  are  considered. 

(2)  Starting  conditions  at  leading  edge  must  be  assumed. 


(3)  The  viscous  condition  of  no  slip  at  the  wing  surface  in  the 
area  of  the  free  vortex  is  not  satisfied  and  this  can  affect 
final  vortex  position  and  strength. 
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(4)  No  secondary  vortex  system  is  included. 

(5)  Vortex  burst  cannot  be  predicted. 


Tools  developed  at  the  Lockheed-Georg ia  Company  have  the  potential  for 
removing  all  of  the  above  limitations.  The  last  three  above  are  viscous 
problems.  To  solve  these,  Lockheed  has  been  under  contract  with  the  Office 
of  Naval  Research  for  approximately  two  years  to  develop  the  "Viscous  Box" 
solution.  References  7 and  8 discuss  previous  efforts  and  the  current 
efforts  are  published  in  this  report.  The  "Viscous  Box"  is  now  a working 
tool.  The  original  efforts  of  the  viscous  box  development  were  directed  to 
a theoretical  investigation  of  spanwise  blowing.  This  concept  controls  the 
leading  edge  vortex  on  straight  as  well  as  swept  wings;  therefore,  much  of 
the  original  analysis  of  References  7 and  8 contains  spanwise  blowing  with 
straight  wings.  This  allowed  the  use  of  a Green's  function  integral  to 
supply  some  of  the  boundary  conditions  and  the  flow-field  velocities.  This 
tool  has  been  retained  to  make  various  parametric  investigations  presented 
in  this  report.  Application  of  the  viscous  box  method  to  delta  and  other 
highly  swept  configurations  forces  the  use  of  an  elliptic  potential  flow 
solution  to  obtain  boundary  conditions.  This  report  discusses  the  integra- 
tion of  the  viscous  box  and  a vortex  lattice  method  as  applied  to  a thin 
delta  wing.  Thick  wings  will  be  treated  in  a subsequent  report. 
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2.  THE  VISCOUS  MODEL 


2.1  Further  Developments  of  the  Viscous  Model 

One  of  the  major  concerns  with  the  numerical  model  for  viscous  flow,  as 
developed  in  Reference  8,  has  been  the  loss  of'  accuracy  associated  with 
the  integration  method.  The  governing  equation  is  that  of  vorticity 
transport, 

V • gradiB  = u>  • grad  V + vV2iB  (1) 

where  u is  the  vorticity  vector  and  is  regarded  as  the  primary  dependent 
variable,  with  V determined  retrospectively  from 

V2V  = - curlSi  . (2) 

In  order  to  maintain  convective  stability  in  the  numerical  integration  scheme, 
it  was  necessary  from  the  outset  to  use  locally  upwind  differencing  for  the 
first  derivative  of  to  in  Equation  (1).  This  resulted  in  first-order  accuracy 
for  integration  in  the  cross-flow  plane.  Such  accuracy  is  acceptable  in  the 
streaming  direction,  since  the  theory  neglects  second  streamwise  derivatives. 
For  the  cross-flow  plane,  however,  this  approximation  can  lead  to  significant 
errors  even  for  moderate  Reynolds  numbers.  In  fact,  the  error,  which  is 
proportional  to  the  product  of  mesh  width  and  local  velocity,  results  in  a 
term  additive  to  second  derivatives  and  thus  is  often  called  "artificial 
viscosity."  In  order  to  upgrade  accuracy  in  the  cross-flow  plane,  first 
derivatives  of  vorticity  with  respect  to  y and  z were  formulated  as  second- 
order  accurate  upwind  differences,  requiring  the  use  of  two  points  upwind  of 
the  current  point.  The  form  of  equation  (l)  upon  expansion  is, 

Fx  = AF  + BF2  + CFy  + DFZZ  + DFyy  (3) 

where  F*(£,n,C),  the  vorticity  vector,  and  A through  D are  3x3  matrices 
involving  the  velocities  and  their  gradients.  Thus,  the  first  y-derivative 
of  the  vector  F is  given  by 


The  basic  integration  molecule  in  use  in  the  present  integration  scheme  is 
depicted  by  the  solid  lines  in  Figure  1.  This  molecule  results  in  a tri- 
diagonal matrix  when  forming  the  implicit  solution  at  each  cross-flow  plane. 
Addition  of  another  upwind  point  results,  In  general,  in  a pent-diagonal 
matrix  for  which  the  inversion  time  is  greatly  increased.  To  avoid  this 
penalty,  it  was  observed  that  since  the  basic  method  is  a fully  nonlinear 
Alternating  Direction  Implicit  (ADI)  scheme,  the  resulting  solution  (assuming 
a convergent  iteration)  should  be  a valid  one  if  values  from  the  previous 
iteration  were  used  for  the  second  upwind  points  (i.e.  at  m±2,  n ± 2) . These 
are  represented  by  the  dashed  lines  in  Figure  1 and  complete  the  augmented 
molecule  for  second-order  accuracy.  This  approach  provides  second-order 
accuracy  with  the  minimum  penalty  in  computing  time. 

The  effect  of  more  accurate  integration  can  be  seen  in  relative  terms 
by  comparing  results  for  a typical  flow  condition.  Figures  2 and  3 
present  a comparison  of  the  ^-component  of  vorticity  and  the  W-component  of 
velocity,  respectively,  as  obtained  by  the  two  methods.  The  flow  condition 
is,  as  usual,  that  of  a vortex  over  a plate  with  cross-flow.  The  vortex  is 
located  between  z = .l»  and  z = .6  and  centered  at  a height  y of  roughly  .**. 

The  effect  of  increased  flow  reversal  is  reflected  in  both  figures.  This 
is  to  be  expected  as  a result  of  removing  artificial  viscosity.  Some 
estimates  of  accuracy  will  be  discussed  in  a later  section. 


Figure  3 Comparison  of  First  and  Second  Order  Accurate  Differencing . W Velocity 
Component  at  X = 3,  A Y = A Z = . 1,  11x21. 


The  Green's  function  method  for  approximating  the  potential  flow  field 
(Reference  8)  has  been  extended  to  represent  plane-of-symmetry  conditions  in 
the  cross-flow  plane.  It  was  concluded  that  some  device  was  required  for 
approximating  potential  flow  conditions  peculiar  to  delta  wings.  Thus, 
highly  swept  delta  wings  exhibit  strong  effects  /near  the  plane  of  symmetry 
due  to  the  vortex  on  the  reflected  leading  edge.  In  earlier  studies,  the 
leading  edge  of  the  surface  was  of  primary  concern  in  modeling  potential 
flow  effects,  and  the  other  edge  of  the  surface  was  envisioned  as  a "down- 
stream" edge.  A more  accurate  representation  for  delta  wings  is  to  consider 
the  y-z  plane  to  be  the  cross-flow  plane  including  the  reflected  leading 
edge.  The  accuracy  of  this  approximation  improves  with  increasing  sweep 
angle.  Figure  b depicts  this  situation,  where  the  computational  domain  is 
now  to  the  right  of  the  centerline  and  the  prescribed  downstream  portion  of 
the  flow  field  is  just  a reflection  about  the  centerline  of  the  flow  to  the 
right.  This  representation  is  conceptually  similar  to  the  classical  conical 
flow  treatment  of  delta  wings  in  potential  flow.  An  attempt  was  made  to 
improve  the  approximation  by  performing  a rotation  of  the  reflected  quantities 
as  depicted  in  Figure  A (b) . However,  when  this  was  applied  to  the  box,  the 
numerical  integration  became  erratic  and  gave  anomalous  results.  This  may  be 
related  to  the  ellipticity  implied  by  the  procedure;  in  any  case,  the  rotation 
was  abandoned  in  favor  of  the  simpler  representat ion.  Figure  5 presents 
typical  results  obtained  with  this  method  for  a 65°  swept  delta  wing  with 
e “1*3  and  a *22°.  Depicted  are  the  ^-component  vorticity  contours  with 
V-W  vectors  superimposed.  A weak  starting  vortex  was  used  and  the  contours 
shown  are  augmented  by  vorticity  convected  in  from  the  leading  edge.  This 
result  demonstrates  the  importance  of  plane-of-symmetry  conditions.  Since 
there  is  no  w-component  of  velocity  on  the  right-hand  side  of  the  domain,  the 
vorticity  fed  in  from  the  leading  edge  is  "trapped"  and  must  accumulate  with 
distance  downstream,  thus  augmenting  vortex  strength. 

2.2  Significance  of  the  Leading  Edge  Vorticity 

The  Green's  function  method  has  provided  a useful  device  for  various 
parametric  studies.  While  the  method  is  not  intended  as  a part  of  the  final 
computational  package,  it  has  aided  in  gaining  insight  regarding  sensitivity 
to  edge  conditions  and  in  determining  the  relationship  between  boundary  layer 
and  vortex.  The  details  of  this  method  are  given  in  Reference  (8). 
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Figure  5 . Velocity  Vectors  and  Contours  of  5 Component  of 
Vorticity  with  Simulated  Plane  of  Symmetry 


Some  computations  were  performed  to  demonstrate  that,  using  the  Green's 
function  method,  the  necessary  leading  edge  vorticity  could  be  obtained  as  a 
part  of  the  boundary  value  solution.  It  appears,  however,  that  this  is  not 
appropriate.  The  problem  of  immediate  concern  is  that  of  a sharp  leading 
edge  surface.  The  discussion  following  should,  ^rowever,  be  valid  whether 
the  leading  edge  is  sharp  or  rounded,  as  long  as  separation  is  very  near  the 
leading  edge.  Examination  of  experimental  results  for  the  V-component  very 
near  the  leading  edge  has  shown  that  the  original  assumption  used  for  V in 
the  Green's  function  was  in  fact  quite  accurate.  The  situation  is  depicted 
in  Figure  6.  The  boundary  layer  detachment  point  is  shown  at  the  forward 
most  point  of  the  leading  edge;  and  in  ’’is  case,  the  boundary  layer  profile 
is  coincident  in  direction  with  V as  <!■..  red  from  the  coordinate  system  of 
the  viscous  box.  However,  within  the  b'  udary  layer  thickness,  6,  the 
actual  profile  shape  can  be  expected  to  differ  from  the  linear  one  assumed 
for  V.  But  the  important  thing  is  to  require  that  the  vorticity  contained 
in  the  boundary  layer,  and  thus  transported  into  the  viscous  box,  be  repre- 
sented on  the  box  boundary.  The  obvious  suggestion  then,  is  to  adjust  the 
thickness  6 assumed  for  the  V-component  so  that  the  vorticity  computed  from 
the  straight  line  segment  of  V is  equal  to  the  average  value  in  the  actual 
boundary  layer.  This  reduces  the  behavior  of  V near  the  leading  edge,  and 
its  associated  boundary  vorticity,  to  the  specification  of  a single  parameter: 
6. 


As  noted  before,  the  apparent  value  of  vorticity  appropriate  to  the 
central  test  condition  of  the  65°  swept  model  is  approximately  A3-  With  the 
computer  program  modified  to  automatically  fix  both  V and  leading  edge 
boundary  vorticity  upon  specification  of  6,  a series  of  runs  were  made  vary- 
ing 5 or  5i.e.  to  create  a range  of  flow  conditions  as  depicted  in  Figure  7- 

(Note:  5].e.  - 3v/8z  across  the  boundary  layer  ahead  of  the  box.)  The 
values  of  e were  20,  30,  A3,  50,  and  60.  For  each,  the  characteristics 
of  the  flow  are  computed  at  the  inlet  face  of  the  viscous  box  and  at  X»1. 

To  illustrate  the  effect  of  this  input,  ^-vorticity  contours  for  X * 1 are 
presented  for  each  of  the  chosen  E}.e.  values  in  Figure  7. 
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VORTICITY  EQUAL  TO  THE 
AVERAGE  VALUE  IN  THE 


Figure  6 . Leading  Edge  Conditions 


In  each  of  the  figures  the  total  £-vorticity,  Z£,  for  the  entire  plane 
at  X*1.0  is  presented  in  the  right-hand  corner  of  the  figure.  There  is  a 
strong  growth  in  I£  as  £i.e.  increases,  and  appears  to  peak  between  £].e.  = 50 
and  60.  The  value  of  real  interest  is  the  sum  of  all  vorticity  off  the  wing 
surface  not  including  £j.e>  . This  has  not  bee^determined  for  these  figures 
but  is  easily  accomplished.  It  is  this  total  and  its  centroid  that  must  be 
used  to  represent  the  free  vortex  system  used  in  the  vortex  lattice  program 
(Section  3*M. 


Results  of  this  investigation  have  shown  that  e must  be  on  the  order 
of  40  to  60.  This  agrees  with  experiments  that  have  shown  8(V/Woo)/3Z  on  the 
order  of  35  to  50.  This  leaves  no  doubt  of  the  importance  of  determining  an 
accurate  leading  edge  boundary  layer  at  the  point  of  separation. 

2.3  Accuracy,  Mesh  Size,  and  Viscosity 

Attempts  have  been  made  to  quantify  the  accuracy  of  the  viscous  box 
method.  This  is  a difficult  task  because  exact  solutions  are  not  available 
for  comparison.  Also,  it  is  not  clear  what,  if  any,  single  parameter  can  be 
used  as  a measure  of  relative  accuracy  between  two  computations  at  different 
mesh  sizes.  It  was  thought  that  an  integral  measure  of  the  primary  vorticity 
component  (i.e.,  £,  the  one  associated  with  the  vortex  strength)  would  suffice 
in  this  regard.  Thus,  as  mesh  size  is  decreased,  an  asymptotic  behavior 
should  become  obvious.  A series  of  three  runs  were  made  in  which  mesh  size 
was  successively  decreased  and  the  integral  measure,  indicated  as  £2  in  Figure 
8,  was  computed  after  one  x-step.  It  is  not  clear  from  the  figure  that  con- 
vergence has  occurred.  Another  measure  was  tried  for  the  same  data:  The 
maximum  value  of  £-vorticity  was  obtained  for  each  run  and  this  is  shown  as 
£max  'n  figure  8.  Here  the  results  are  more  encouraging.  As  expected,  the 

maximum  vorticity  increased  with  decreasing  mesh  size  and  it  appears  to  be 

\ • 

asymptotic.  The  slightly  greater  value  of  Smax  at  the  intermediate  mesh 
size  may  be  due  to  an  inconsistent  distribution  of  initial  value  vorticity 
between  the  intermediate  (1/Ay  = 15)  and  final  (1/Ay  = 20)  mesh  sizes. 

Another  test  of  accuracy  lies  in  determining  the  lower  limit  of  viscosity 
coefficient  for  which  no  further  change  in  flow  field  properties  can  be 
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Figure  8 Effect  of  Mesh  Size  on  0 « JJ  Sdydz,  and 
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observed.  A series  of  three  runs  were  made  varying  only  the  viscosity  coeffi- 

—2  -3  -4 

cient.  Values  of  v ■ 10  ,10  , and  10  were  used.  Again,  the  quantity, 

previously  referred  to,  was  computed  and  this  is  shown  in  Figure  9-  This 
figure  indicates  that  the  asymptotic  value  of  v lies  between  10“3  and  10“4 
when  mesh  size  is  11  x21.  The  implication  of  tl^s  is  that  at  the  mesh  size 
which  is  typically  used  in  our  computations,  the  upper  limit  of  Reynolds 
number  dictated  by  v lies  somewhere  between  10-3  and  10"4.  Beyond  this  point 
the  real  viscosity  is  swamped  by  numerical  inaccuracy.  Of  course,  this 
Reynolds  number  limit  can  always  be  increased  by  reducing  mesh  size.  A 
further  measure  of  the  effect  of  v is  indicated  in  Figure  10.  This  is  a 
comparison  of  a set  of  W-profiles  at  a particular  point  in  the  flow  (z  = .5) 
for  a change  in  viscosity.  The  major  difference  occurs  at  the  outer  edge  of 
the  box  (the  Green's  function  technique  for  velocities  was  used  in  these 
runs)  and  is  seen,  as  expected,  to  be  less  between  the  two  smaller  values  of 
v. 


A spatially  variable  viscosity  coefficient  was  tested  in  the  computer 
code,  and  the  numerical  results  were  found  to  be  stable.  However,  the 
gradient  terms  due  to  variable  viscosity  were  not  included.  If  in  the  future 
a turbulence  model  is  incorporated  in  the  code,  then  it  will  be  necessary  to 
add  these  terms.  The  present  investigation  was  intended  only  to  demonstrate 
feasibility,  and  no  use  has  been  made  of  the  variable  viscosity  feature. 

2. A Effect  of  Sweep  and  Flow  Angularity  at  the  Leading  Edge 

Parametric  studies  were  performed  to  evaluate  the  effects  of  sweep  and 
flow  angularity  on  leading  edge  vortex  behavior.  The  effects  of  both  are 
very  similar;  the  lower  limit  of  sweep  angle  for  which  a vortex  can  be  sus- 
tained is  somewhat  larger  than  25°.  The  flow  is  definitely  unstable  at  that 
angle  and  is  probably  not  stable  at  considerably  higher  angles  if  the  calcu- 
lation is  marched  far  enough  outboard.  The  parametric  study  was  performed 
using  only  one  x-step. 

The  effect  of  flow  angularity  is  a more  interesting  case  and  will  be 
discussed  in  some  detail.  This  effect  is  best  discussed  in  conjunction  with 
Figure  11.  The  vorticity  feeding  the  leading  edge  of  the  vorticity  box  is, 
as  has  been  discussed  in  the  past,  associated  with  the  separating  boundary 
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Figure  10  Effect  of  Viscosity  Coefficient  on  a Typical  Velocity  Profile 
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region  just  upstream  of  the  box.  The  development  of  this  region  is  very  much 
dependent  on  leading  edge  shape,  and  its  effect  on  vortex  development  within 
the  box  is  a function  not  only  of  the  magnitude  of  vorticity  emanating  from 
the  region  but  also  its  direction.  Thus,  when  the  vorticity  vector  is 
parallel  to  the  leading  edge,  as  shown  in  Figure^ll,  a sharp  leading  edge  is 
implied.  As  thickness  is  added  to  the  leading  edge,  it  is  reasonable  to 
expect  that  some  downstream  turning  of  the  boundary  region  will  be  reflected 
through  the  rotation  of  the  surface  streamline.  This  also  implies  a rotation 
of  the  vorticity  which  feeds  the  box.  To  estimate  the  effect  of  such  a 
rotation  on  leading  edge  vortex  behavior,  several  numerical  experiments  were 
performed  in  which  the  magnitude  of  vorticity  was  held  constant,  but  was 
apportioned  between  £ and  e;  components  at  the  leading  edge  so  as  to  reflect 


the  vector  rotation. 


Again,  as  in  the  case  of  sweep,  the  flow  progresses  toward  unstable  con- 
ditions as  the  vector  is  rotated.  Figure  12  summarizes  the  behavior  for  a 
rotation  of  40°  away  from  the  leading  edge  for  a sweep  angle  of  65°.  The  U- 
component  of  velocity  shown  is  for  the  most  critical  z-station,  i.e.  where 
the  first  flow  approaches  reversal.  After  three  x-steps  the  computations 
diverge  due  to  reversal.  This  presumably  corresponds  to  a flow  condition 


for  which  a vortex  cannot  be  sustained. 


2.5  Computation  of  Pressure  in  the  Viscous  Region 


It  is  well  known  in  computational  fluid  mechanics  that  computing  pressure 
in  a viscous  flow  is  very  sensitive  to  numerical  differencing  errors.  The 
procedure  used  here  is  as  follows:  The  momentum  equations  are 


grad  + vV2v  * v 
P 


grad  (^j-) 


V2v  = - curl  to 


then, 


grad  (£■  + j v2)  = -v  curl  <L  + v x 5 . 


Let 
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H(x,y,z)  -J  + lv2 
P 2 


so  that, 

p ■=  p(H  - j v2)  and  Cp  « 

OD 


Then  any  one  of  the  three  equations, 
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may  be  used  to  obtain  Cp.  The  second  of  these  was  chosen  for  integration  in 
the  direction  normal  to  the  surface,  beginning  on  the  outer  edge  of  the  box 
where,  if  the  box  is  large  enough,  there  exists  H*H  . The  equation  is  then 
integrated,  to  second-order  accuracy,  through  the  box  to  the  wall.  The 
pressure  is  computed  at  each  mesh  point;  Figure  13  is  a typical  example  of 
pressure  coefficient  contours  obtained  in  this  manner.  Pressure  at  the 
surface  is  the  important  quantity  and  in  order  to  compute  this  it  is  neces- 
sary to  march  through  all  the  mesh  points  from  above,  thus  accumulating 
numerical  error  at  the  wall.  However,  using  either  of  the  other  two 
equations  leaves  some  doubt  about  the  correct  edge  value  for  H.  Evaluation 
of  accuracy  by  the  present  method  and  several  alternative  methods  is  a 
subject  of  continuing  study.  Detailed  calculations  for  wall  pressure  will 
be  discussed  in  Section  3.8.3* 

2.6  Iteration  and  Step  Length  Sensitivity 

The  major  difficulty  with  the  parabolized  flow  approximation  as  used  in 
the  viscous  box  is  sensitivity  to  initial  conditions.  This  is  expected  of 
parabolic  systems  when  the  marching  variable  is  not  far  removed  from  its 
starting  value.  When  this  variable  becomes  large,  however,  the  solution 
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gure  13  Pressure  Coefficient  Contours 


"forget"  the  initial  conditions  and  become  dominated  by  the  boundary  values. 

If  the  initial  conditions  are  not  physically  correct,  it  becomes  very  diffi- 
cult and  sometimes  impossible  to  converge  the  iteration  for  the  first  step 
away  from  the  starting  plane.  In  the  present  study  of  leading  edge  vortex 
flows,  exact  starting  values  for  some  pos i t ions  .along  the  leading  edge  have 
not  been  available.  The  procedure  in  all  the  computations  has  been  to  make 
a guess  of  starting  conditions  and  step  a distance  generally  of  the  order  of 
box  dimensions  in  the  cross-flow  plane.  At  this  position  a large  number  of 
iterations  are  generally  required  to  converge  to  an  acceptable  three- 
dimensional  solution.  Under  these  conditions,  the  iteration  is  not  conver- 
gent when  the  initial  x-step  is  reduced  to  one  order  of  magnitude  lower  than 
the  lateral  dimensions.  This  simply  means  that  the  numerical  solution 
procedure  is  more  sensitive  to  the  inaccuracy  of  the  assumed  starting  values. 

It  is  noted,  however,  that  once  a valid  solution  is  established,  the  marching 
procedure  is  well  behaved,  both  with  respect  to  stepping  distance  and  to 
change  of  boundary  conditions  with  x. 

A "distance  relaxation"  computation  was  performed  to  determine  if  this 
would  produce  an  improved  starting  solution.  The  analogy  here  is  with  the 
time  relaxation  of  the  full  parabolic  Navier-Stokes  equations,  in  which  steady 
flow  solutions  to  the  (elliptic)  system  are  obtained  by  advancing  the  parabolic 
solution  over  a long  time.  In  the  present  case,  a single  large  x-step  can  be 
used  to  relax  the  solution  independent  of  x.  Since  the  solution  procedure  is 
implicit  and  fully  nonlinear  (i.e.  no  linearization  about  coefficients  at  the 
previous  step),  only  one  step  is  required.  The  conditions  for  which  this  was 
performed  were  those  applying  to  the  first  box  on  the  delta  wing  demonstration 
case,  discussed  in  another  section.  The  point  here  is  simply  to  demonstrate 
the  results  of  such  a computation.  The  box  dimensions  were  of  order  unity 
and  a step  length  Ax  = 106  was  used.  A typical  iteration  history  of  vorticity 
appears  in  Figure  1A.  As  shown,  the  solution  required  16  iterations  to  con- 
verge and  the  iteration  history  is  well  behaved.  In  order  to  perform  the 
computation,  a set  of  boundary  values  were  required  as  well  as  initial 
conditions.  The  boundary  values  of  velocity  were  those  obtained  from  vortex 
lattice  calculations  for  the  starting  face  of  the  first  box,  and  the  leading 
edge  vorticity  was  that  estimated  for  the  first  box  (as  discussed  in  a later 
section).  The  cross-flow  velocity  vector  and  ^-component  vorticity  fields 
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resulting  from  relaxation  with  those  imposed  values  are  shown  in  Figure  15. 
These  would  appear  to  represent  reasonable  initial  values  for  beginning  the 
downstream  march.  However,  the  U-component  of  velocity,  as  seen  from  Figure 
16  has  attempted  reversal.  Actual  reversal  is  not  allowed  in  the  computer 
code  since  this  leads  to  a divergent  iteration.  .-The  profiles  at  z = .15  and 
.25  are  representative  of  the  behavior  just  forward  of  the  apparent  vortex 
center,  and  without  constraint  these  profiles  would  indicate  reversal.  The 
profiles  near  the  vortex  center  again  show  the  characteristic  axial  accelera- 
tion under  the  vortex.  Apparently,  the  given  boundary  conditions  will  not 
allow  a vortex  to  be  sustained  over  a very  long  distance.  The  indicated 
reversal  could  be  associated  with  vortex  bursting,  though  no  axial  pressure 
gradient  is  present.  Further  work  is  needed  to  improve  starting  values. 

This  can  result  in  significant  reduction  of  computer  time  since,  with  an 
accurate  solution  established,  succeeding  steps  require  fewer  iterations. 

The  results  depicted  in  Figure  17  clearly  demonstrate  the  advantage  of 
an  accurate  starting  solution.  This  iteration  history  resulted  for  the  same 
initial  and  boundary  conditions  discussed  above  but  with  an  X-step  of  unity. 
The  solid  lines  show  the  convergence  pattern  on  successive  steps  following 
the  first  one,  which  required  16  iterations  to  converge.  These  later  steps, 
which  are  begun  using  the  solution  from  the  previous  step  as  the  initial 
guess  for  beginning  the  iteration,  are  seen  to  converge  in  about  three 
iterations.  In  general,  a few  more  iterations  than  this  would  be  expected 
because  of  changing  boundary  conditions  at  each  X.  Figure  17  also  reveals  a 
difficulty  with  the  iteration  procedure  resulting  when  the  vorticity  is  near 
zero.  Although  the  "shock"  occurring  in  the  iteration  does  not  destabilize 
the  process,  it  may  well  delay  convergence.  The  source  of  this  behavior  has 
not  been  isolated  yet,  but  no  difficulty  is  anticipated  in  correcting  it.  The 
dashed  line  in  Figure  17  depicts  the  iteration  history  at  X * 1 , resulting 
from  first-order  accurate  differencing  of  the  wall  value  of  vorticity.  The 
procedure  used  in  the  viscous  box  code  has  always  been  to  update  wall 
vorticity  using  a second-order  accurate  one-sided  first  difference.  It  has 
been  noted  that  such  differencing  does  not  conserve  vorticity  whereas,  at 
least  theoretically  first-order  differencing  does.  Our  numerical  experi- 
ments do  not  provide  sufficient  evidence  to  support  a change  to  the  lower 
accuracy  at  present.  Further  experiments  are  planned  to  examine  this 
question. 


26 


■ 


K,  ITERATION  NUMBER 


Figure  17  Effect  of  Initial  Values  on  Iteration  History,  £ at  y = 0 


3.  INTEGRATION  OF  THE  VISCOUS  BOX  AND 
ELLIPTIC  POTENTIAL  FLOW  MODELS 


3 . 1 General 

The  purpose  of  the  current  efforts  to  integrate  the  potential  and 
viscous  flow  models  is  to  determine  the  problem  areas  involved  and  ascertain 
the  feasibility.  A complete  converging  solution  has  not  been  accomplished 
at  th.is  time,  but  sufficient  effort  has  been  made  to  show  the  feasibility  of 
the  integration  for  thick  wings  with  rounded  leading  edges  as  well  as  sharp 
thin  wings. 

To  test  the  concept,  a delta  wing  was  chosen  as  the  configuration  to 
analyze.  The  concept  is  not  limited  to  delta  wings,  but  much  is  already 
known  about  the  delta  wing  flow  at  least  qualitatively  if  not  quantitatively. 
Some  discussions  of  the  technique  applied  to  wings  with  leading  edge  strake 
is  presented  subsequently. 

The  techniques  presented  below  were  not  without  problems;  however,  these 
problems  proved  to  be  valuable  in  understanding  some  of  the  more  important 
parameters.  These  parameters  along  with  correct  and  incorrect  methods  of 
application  are  discussed. 

Although  a single  converging  computational  program  is  not  yet  available, 
the  tools  that  have  been  developed  can  be  put  together  to  provide  a viscous 
flow  solution  to  configurations  containing  interacting  free  vortices  and 
lifting  surfaces.  These  tools  will  allow  solution  of  high  angle-of-attack 
maneuvering  aircraft  as  well  as  other  vortex  flow-field  problems  such  as  up- 
swept  fuselage  afterbody  drag  and  upper  surface  blown  wings,  both  chordwise 
and  spanwise  blowing. 

3.2  Modeling  Technique 

The  technique  of  solving  the  vortex/wing  interaction  problem  with  a 
hybrid  integration  of  viscous  and  potential  flow  models  is  now  clear.  A flow 
diagram  is  presented  in  this  section.  Two  potential  flow  elliptic  solutions 
must  be  integrated  with  two  viscous  parabolic  solutions  to  form  a hybrid 
solution  containing  the  best  of  each  model.  The  four  models  are: 
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(1)  Thick  or  thin  wing  potential  flow  model 

(2)  Free-vortex  system 

(3)  Viscous  box 

(A)  Leading  edge  boundary  layer. 

Figure  18  shows  a flow  diagram  for  integrating  these  models.  At  this  time 
only  the  thin  wing,  sharp  leading  edge  delta  wing  has  been  investigated.  In 
general,  however,  the  four  methods  are  combined  as  follows.  After  a suitable 
first  estimate  of  the  free  vortex  system  (strengths  and  positions,  section 
3.3),  the  free-vortex  system  is  combined  with  the  thick  or  thin  wing  potential 
flow  model.  The  model  provides  the  pressure  distribution  around  the  leading 
edge  as  well  as  boundary  velocities  on  the  viscous  box  boundaries.  The 
pressure  distributions  are  used  in  the  leading  edge  boundary  layer  solution 
to  get  vorticity  fed  to  the  viscous  box.  With  initial  boundary  conditions 
defined,  the  viscous  box  is  then  used  to  find  more  exact  locations  and 
strengths  of  the  free-vortex  system.  Then  the  procedure  is  repeated  until 
convergence  occurs.  Convergence  occurs  when  the  total  strength  and  centroid 
of  the  circulation  in  the  viscous  box  do  not  change  more  than  certain  pre- 
scribed limits.  These  procedures  are  detailed  in  Section  3-2.2. 

The  above  procedures  have  been  applied  to  a simple  thin  sharp-edge  delta 
wing  with  65°  leading  edge  sweep  at  22°  angle  of  attack.  This  angle  was 
chosen  because  it  should  have  the  leading  edge  vortex  with  incipient  burst 
near  the  wing  trailing  edge.  Figure  19  shows  a planform  and  projected  view 
of  the  wing  with  viscous  boxes  in  place.  Note  that  the  front  of  the  box 
starts  at  the  leading  edge.  In  a subsequent  section  it  is  noted  the  box 
should  extend  slightly  ahead  of  the  leading  edge  of  the  wing.  More  boxes 
would  be  desirable  for  detailed  solution,  but  to  keep  computer  time  down  while 
developing  the  procedure,  only  six  boxes  are  used.  The  free-vortex  system 
from  the  wing  apex  to  the  first  box  should  be  extrapolated  according  to  the 
relations  used  for  approximating  the  initial  vortex  system  (see  3-2).  This 
assumption  will  not  significantly  affect  the  overall  results  or  purpose  of 
this  effort. 
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Figure  18  Flow  Diagram  for  Wing  Vortex  Viscous/Potential  Flow  Interaction  Models 
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3.2.2  Detailed  Procedures 


The  following  step-by-step  procedure  is  recommended  for  obtaining  a 
solution  to  the  vortex/wing  interaction  problem. 

(1)  Estimate  leading  edge  free-vortex  strength  and  core  position 
from  approximate  method  described  in  Section  3-3. 

(2)  Define  thick  wing  potential  flow  paneling  to  give  good 
pressure  distribution  around  leading  edge.  For  the  thin  wing 
a smaller  number  of  panels  may  be  used,  but  be  sure  not  to 
place  viscous  boxes  near  bound  vortex  elements. 

(3)  Determine  number  and  size  of  viscous  vortex  boxes  to  contain 
all  of  vorticity  of  step  (1)  above.  For  thick  wings,  start 
box  on  an  element  aft  of  leading  edge  to  allow  for  boundary 
layer  growth  and  leading  edge  curvature.  See  Section  3-6. 
Start  box  near  collocation  points  at  leading  edge.  For  thin 
wings,  start  box  ahead  of  leading  edge. 

(A)  Provide  card  or  tape  output  from  potential  flow  model  of 
collocation  points. 

(5)  Set  up  free-vortex  model  with  bound  vortex  elements  compatible 
with  potential  flow  model. 

(6)  Use  free-vortex  model  to  compute  (a)  downwash  at  potential 
flow  model  collocation  points,  (b)  velocities  where  surface 
pressures  are  to  be  calculated,  and  (c)  velocities  at 
specified  edges  and  grid  points  on  the  faces  of  the  boxes. 

Put  output  on  cards  or  tape  to  supply  to  potential  flow  model. 

(7)  Using  the  downwash  calculations  from  step  (6),  solve  the 
potential  flow  model  system  and  compute  (a)  velocities  where 
surface  pressures  are  to  be  calculated  and  (b)  velocities  at 
the  specified  edges  and  grid  points  on  the  faces  of  the  boxes. 


, V-  • . . 
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(8)  Add  velocities  of  steps  (6)  and  (7)  for  surface  pressures  and 
box  locations  to  get  total  three  component  velocities. 

(9)  Convert  surface  velocities  to  pressures  and  solve  the 
boundary  layer  equations  for  flow  around  leading  edge  into 
viscous  boxes. 

(10)  Beginning  with  the  first  viscous  box,  use  the  viscous  box  method 
to  compute  the  flow  field  in  the  box,  at  its  face,  its  lee  end 
and  desired  stations  in  between. 

(11)  Use  vorticities  and  velocities  computed  in  step  (10)  for  the  lee 
end  of  box  1 for  input  to  box  2 where  box  1 covers  box  2.  See 
Section  3-5- 

(12)  Where  the  face  of  box  2 is  not  covered  by  the  lee  end  of  box  1, 
use  the  velocities  computed  from  the  free  vortex  system  and 
potential  flow  model  for  input  to  the  face  of  box  2.  See 
Section  3-5- 

(13)  Repeat  steps  (10),  (11),  and  (12)  for  all  of  the  remaining  boxes. 

(1A)  At  each  station  where  viscous  box  calculations  are  made,  integrate 
the  vorticity  to  obtain  the  total  circulation  about  the  box  and 
its  centroid. 

(15)  Divide  the  box  at  each  station  into  approximately  negative  and 
positive  vorticity  and  find  the  centroid  and  strength  of  circu- 
lation for  each  area.  See  Section  3*7- 

( 1 6)  Define  a new  free-vortex  system  based  on  the  distributed  circu- 
lation computed  in  step  15.  See  Section  3.k. 

/ 

/; 

(17)  Repeat  steps  6 through  14  and  check  to  see  if  total  circulation 
and  centroid  changes  are  within  prescribed  limits  (say  2 or  3%). 


A 


( 1 8 ) If  step  (17)  converges,  calculate  overall  forces  and  moments 
and  pressure  distribution  with  the  combined  free-vortex  model 
of  step  ( 1 6)  and  the  thick  wing  potential  flow  model.  Pressures 
at  collocation  points  should  be  sufficient. 

The  following  sections  provide  additional  detail  of  the  above  steps. 

3.3  Approximation  for  Estimating  Vortex  Strength  and  Position 

A correlation  of  data  from  various  experiments  and  the  theory  of  Kandil, 
et  al.^,  has  been  made  to  provide  the  vortex  position  as  a function  of  the 
parameter  "a"  = a/tany , where  a is  the  angle  of  attack  and  y is  the  angle  made 
by  the  wing  leading  edge  and  the  root  chord.  Figure  20  shows  this  correlation 
of  various  experimental  data^*1^.  Some  of  the  various  values  of  X/Cr  are 
given.  There  is  some  scatter  due  to  X/Cr  varying  but  does  not  seem  consis- 
tent. Therefore,  a line  is  put  through  the  points  with  the  conical  flow 
Smith  theory  as  a guide11.  The  value  of  "a"  for  the  wing  we  have  chosen  to 
analyze  is  0.8234  (a  =22°).  This  is  in  between  the  two  sets  of  data  shown  in 
Figure  20. 

To  determine  the  effect  of  "a",  a crossplot  is  made  in  Figure  21.  Again 
using  Smith  as  a guide,  lines  are  drawn  for  y and  z as  a function  of  "a". 

This  is  crude  for  more  values  of  "a"  are  needed  from  experiment.  For  our 
value  of  "a"  the  curve  of  Figure  20  gives  good  agreement  with  the  vertical 
location  of  the  vortex  found  in  our  laser  velocimeter  tests. 

Vortex  strength  is  determined  as  a function  of  X/Cr,  a and  y and  is 
plotted  in  Figure  22.  Theoretical  values  for  two  delta  wing  configurations, 
analyzed  in  Reference  6,  have  been  correlated  to  give  a universal  value  Kv 
as  a function  of  X/Cr,  where 


U00S  sin2a  cosa 


Note,  Kvsin2a  cosa  = C(_v  as  given  by  Polhamus,  Reference  3,  in  his 
suction  analogy.  This  is  not  the  same  Kv  but  the  form  and  purpose  is  similar. 
It  should  be  noted  that  the  f obtained  from  Figure  22  does  account  for  the 
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Figure  21  Vortex  Location  for  Flat  Delta  Wings 


Figure  22  Delia  Wing  Leading  Edge  Vortex  Strength  Distribution  (Thin  Sharp  Edge  Deltas) 


Kutta  condition  at  the  wing  trailing  edge,  and  this  contrasts  with  Smith's 
conical  flow  values.  Figure  23  shows  the  values  of  r computed  for  the  delta 
wing  used  for  the  investigation  of  this  report.  Use  of  these  values  are 
discussed  in  more  detail  in  the  following  section. 

3.3-1  Suggested  Technique  for  Estimating  Free  Vortex  Strengths 
and  Positions  on  the  Wing  With  Leading  Edge  Strake 

This  effort  is  suggested  to  illustrate  the  application  of  the  hybrid 
method  to  a very  complex  but  apparently  solvable  problem.  Figure  24  shows 
the  planform  of  a semi-span  wing  recently  tested  at  low  subsonic  speeds. 
Vortex  position  and  strength  for  the  outer  basic  wing  are  approximated  in 
the  same  manner  as  the  delta  wing.  An  equivalent  delta  is  outlined  as  dotted 
lines  for  both  the  outer  wing  and  the  strake.  It  must  be  remembered  that 
this  is  only  the  first  approximation  and  the  actual  vortex  position  and 
strength  will  eventually  be  calculated  with  the  hybrid  method.  It  is  of 
interest  to  note  that  the  approximate  vortex  position  as  shown  in  Figure  24 
appears  very  close  to  surface  flow  observations. 

The  trailing  strake  vortex  aft  of  the  strake  will  not  increase  in 
strength  over  the  rest  of  the  wing.  The  position  of  this  vortex  element  will 
have  to  be  determined  by  tracking  methods  similar  to  those  used  by  Kandil  or 
Rao  in  References  5 and  6. 

A very  limited  attempt  was  made  under  this  contract  to  simulate  the 
strake  vortex  effect  on  the  outboard  vortex  using  the  viscous  box  as  shown  in 
Figure  24.  A preliminary  computation  was  performed  to  evaluate  the  effect 
of  a strake  vortex  on  a vortex  formed  over  a wing  of  moderate  sweep  angle.  A 
box  calculation  was  performed,  which  was  initialized  with  a vortex  plus  a 
velocity  field  due  to  the  "upstream"  strake  vortex.  The  downstream  boundary 
conditions  were  those  appropriate  to  a Green's  function  approximation  of  a 
surface  of  25°  sweep  and  22°  angle  of  attack  plus  edge  velocities  due  to  the 
strake  vortex.  The  vortex  was  immediately  washed  out,  after  one  step  in  the 
spanwise  direction.  The  apparent  conclusion  to  be  drawn  is  that  a fully 
interacting  lifting  surface  calculation  would  be  required  in  order  to 
generate  the  necessary  velocity  field  on  the  edges  of  the  box.  Thus,  the 
lattice  method  would  be  applied  to  the  strake-plus-wing  with  two  vortices 
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Figure  2*»  Wing  with  Stroke  and  Free  Vortices 


imposed  at  approximately  the  correct  locations  and  with  prescribed  strength 
distributions.  With  this  information,  plus  information  about  the  separating 
region  upstream,  the  box  computation  could  be  expected  to  estimate  the 
viscous  development  in  the  vortex  formation  region. 

3.^  Vortex  Lattice  Method 

For  the  thin  delta  wing  investigated  in  the  feasibility  study,  the  flow 
has  been  modeled  with  a vortex  lattice  representation  of  the  wing  surface  and 
with  two  alternative  vortex  system  representations  of  the  separated  leading 
edge  vortex.  These  provide  initial  conditions  for  input  into  the  vorticity 
box  program,  which  performs  a viscous  calculation  with  the  Navier  Stokes 
equations  over  rectangular  box  faces  normal  to  the  wing  surface,  see  Figure 
19*  First,  the  velocities  induced  by  the  leading  edge  vortex  system  are 
calculated  at  the  collocation  points  on  the  wing  and  on  the  faces  and  edges 
of  the  boxes.  Then,  the  velocities  induced  by  the  wing,  under  the  influence 
of  the  leading  edge  vortex,  are  found  for  the  faces  and  edges  of  the  boxes 
and  superimposed  on  those  induced  by  the  leading  edge  vortex  system.  This 
provides  a starting  condition  from  which  the  vortex  strength  and  locations 
can  be  modified  via  the  output  of  the  vorticity  box  program.  Output  of  the 
vorticity  box  is  discussed  in  Section  3.7. 

The  vortex  lattice  arrangement  is  shown  in  Figure  25;  it  features  a 
spanwise  cosine  spacing  of  nine  panels  and  five  equally  spaced  chordwise 
panels.  The  spanwise  panel  spacing  concentrates  them  at  the  leading  edge 
and  near  the  inboard  edge  of  the  rolled-up  vortex  sheet. 

The  first  of  the  leading  edge  vortex  systems  is  shown  in  Figure  26. 

There  are  five  vortex  systems  in  all  that  progressively  feed  into  the  core 
vortex,  the  location  and  strength  of  which  was  defined  using  a conical 
assumption  and  guidance  from  experimental  data  (Section  3-3).  The  individual 
vortex  systems  consist  of  vortex  segments  which  are  input  in  order  from  the 
root  to  the  leading  edge,  from  the  leading  edge  to  the  core,  and  then  to  the 
core  positon  at  the  trailing  edge  and  thence  downstream.  The  root  to  leading 
edge  segment  is  aligned  with  the  wing  bound  vortex. 


Results  using  this  system  are  shown  in  Figures  27  and  28.  The  spanwise, 
vertical,  and  axial  velocity  perturbations  are  presented  for  the  top  edge  of 
the  second  box.  Most  remarkable  is  the  axial  induction  of  flow  for  the 
leading  edge  vortex  system.  Also,  velocity  components  are  calculated  for  the 
inboard  and  outboard  vertical  edges  of  the  box  and  demonstrate  a strong  upwash 
and  strong  lateral  velocities  on  the  outboard  edge. 

A second  and  more  realistic  leading  edge  vortex  system  which  distributes 
vorticity  in  the  feeding  sheet  is  shown  in  Figure  29.  There  are  two  vortex 
systems  at  each  chordwise  position,  one  of  which  feeds  into  the  core  and  the 
other  of  which  feeds  into  the  vortex  sheet.  The  filaments  feeding  into  the 
sheet  are  arbitrarily  assumed  to  feed  into  the  vortex  core,  and  the  vortex 
systems  at  each  chordwise  position  has  the  same  total  strength  as  that  assumed 
for  the  concentrated  vortex  system.  Also,  the  vortices  are  input  in  a manner 
similar  to  that  for  the  concentrated  vortex,  except  that  more  vortex  segments 
must  be  used  to  model  the  sheet.  For  the  first  pass,  the  filaments  compris- 
ing the  vortex  sheet  are  assumed  to  be  distributed  along  a radial  arc  from  the 
leading  edge  to  a position  directly  over  the  vortex  core.  This  sheet  position 
is  determined  later  by  the  vorticity  box.  Once  the  vortex  sheet  filaments 
leave  the  chord  position  at  which  they  are  introduced,  they  are  assumed  to 
follow  a straight  conical  line  to  their  position  on  the  trailing  edge  vortex 
sheet. 


Figures  27  and  28  show  the  velocity  components  for  this  leading  edge 
vortex  system  and  are  qualitatively  similar  to  those  of  the  concentrated 
vortex  system.  The  most  remarkable  difference  is  the  stronger  vertical  compo- 
nent induced  on  the  outboard  edge  of  the  box  due  to  the  proximity  of  the 
vortex  sheet  filaments.  Experimental  data  on  these  plots  are  discussed  in 
Section  4.3. 


3 . 5 Interface  Between  Different  Size  Vorticity  Boxes 

In  order  to  apply  the  vorticity  box  calculation  over  a delta  wing  plan- 
form,  it  is  necessary  to  vary  box  size  while  marching  in  the  downstream 
direction,  as  indicated  in  Figure  19*  Since  the  vortex  grows  with  distance 
downstream,  the  physical  dimensions  of  the  box  must  increase  in  order  to 
contain  the  motion  free  of  boundary  interference.  The  need  for  a cirvi linear 


coordinate  system  to  treat  this  configuration  is  obvious.  Such  a formulation 
will  be  discussed  in  a later  section.  The  present  objective  is  to  test  the 
feasibility  of  coupling  the  viscous  box  with  a potential  flow  model,  and  for 
this  purpose  the  rectangular  domain  will  suffice.  The  major  requirements  for 
implementing  the  coupling  procedure  are  to  read  in  boundary  values  from  the 
lattice  method  and  to  perform,  at  each  box  interface,  an  interpolation  of 
velocities  and  vorticities  to  the  new  grid  points.  The  vortex  lattice 
program  generates  velocities  in  a coordinate  system  as  shown  in  Figure  19* 
These  velocity  components  must  be  related  through  the  sweep  angle  and  angle 
of  attack  in  order  to  conform  to  the  box  coordinate  system.  The  rotations 
and  appropriate  redefinitions  of  velocity  components  are  as  follows: 


(15) 


where  (freestream  direction),  (spanwise),  and  (vertical)  are  the 
velocity  components  as  defined  in  the  lattice  program.  Redistribution  of  the 
solution  to  a new  grid  is  accomplished  by  interpolating  along  vertical  lines 
and  interpolating  a second  time  along  horizontal  lines.  The  grid  lines  are 
depicted  in  Figure  30.  Thus,  if  Q.  represents  any  one  of  the  six  components 
to  be  interpolated,  the  procedure  is  as  follows: 


cosa  sinA 

cosA 

-sina  sinA 

Uj i’ 

= «i 

s i na 

0 

cosa 

p * 

V* 

[cosa  sinA 

-sinA 

-sina  cosA, 

Qni, j = Qmj  + 


fQj>j  - Qi-ul 


Vi  ■ Vi-1 


(ym  - y i-i ) » ■ = i»  m 


(16) 


where  yj_^  <ym<yi*  Indices  i and  j refer  to  the  old  mesh  and  m and  n refer 
to  the  new  mesh.  Equation  (l6)  interpolates  along  vertical  lines;  the  inter- 
polation along  horizontal  lines  then  gives  the  final  result: 


Vn  “ Qfn.j-I  + ( Zj  - ‘ ^Zn  " zj-1^»  ” 1 ’ 


(17) 


where  Zj_  < zn  < z j . Thus,  as  is  clear  from  Figure  30,  the  upstream  solution 
is  distributed  to  the  new  box,  but  its  dimensional  extent  is  unchanged.  This 
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leaves  a portion  of  the  new  box  without  initial  values.  Velocities  for  this 
region  are  taken  from  the  lattice  program.  The  vorticity  for  this  region  is 
left  zero  except  at  the  wall,  where  the  no-slip  condition  must  be  satisfied. 

A starting  value  for  the  wall  is  computed  from  the  given  velocity  field  with 
velocity  at  the  wall  assumed  to  be  zero. 

A better  arrangement  in  shifting  to  a larger  box  would  be  to  maintain  a 
constant  mesh  and  simply  increase  the  number  of  mesh  points  in  the  larger 
box.  This  would  improve  the  accuracy  of  the  numerical  integration  and  avoid 
the  interpolation  requirement.  However,  the  computing  cost  becomes  prohibi- 
tive in  this  case,  thus  the  only  viable  alternative  presently  is  to  hold 
constant  the  number  of  mesh  points.  For  purposes  of  the  present  demonstra- 
tion of  mesh  was  held  at  N = 11,  M = 21  , so  that  the  depth  of  the  box  is  always 
one-half  the  width. 

The  procedure  for  computing  one  pass  over  the  wing  surface  can  be 
briefly  summarized  in  a flow  chart  as  shown  in  Figure  31.  The  upstream  face 
of  the  first  box  is  initialized  with  estimated  or  measured  values  of  velocity 
and  vorticity.  The  solution  is  then  marched  downstream.  If  there  is  more 
than  one  x-step  in  a particular  box,  the  appropriate  edge  values  of  velocity 
are  read  in  from  cards  (or  from  some  other  storage  device)  at  each  y-z  plane 
and  the  boundary  values  of  vorticity  are  similarly  read  in.  For  each 
succeeding  box  face,  both  interior  and  boundary  values  of  velocity  are  read 
from  cards  and  an  overwrite  is  performed  by  the  interpolation  procedure  for 
part  of  the  new  box  face.  The  solution  proceeds  in  this  manner  to  the  output 
of  the  final  box.  At  each  x-step  the  solution  is  written  to  a storage  device, 
to  be  accessed  later  for  further  processing  such  as  automatic  plotting  and 
preparation  of  the  next  lattice  iteration. 

3.6  The  Leading  Edge  Problem 

The  leading  edge  vortex  strength  and  position  has  been  found  to  be 
strongly  affected  by  the  nature  of  the  leading  edge  separation.  Figure  32 
shows  a typical  effect  of  changing  the  leading  edge  radius  on  a 67.1°  swept 
delta  wing1^.  It  is  noted  that  the  increasing  leading  edge  radius  causes  a 
large  increase  in  chordwise  thrust,  Cs,  with  only  a small  loss  in  normal 
force,  Cfj,  or  vortex  lift.  Also,  the  aerodynamic  center  is  moving  aft 
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Figure  31  Procedure  for  Advancing  the  Viscous  Solution 
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indicating  a shift  aft  of  the  leading  edge  vortex  with  increasing  leading 
edge  radius.  This  movement  of  the  vortex  core  is  associated  with  the  orienta- 
tion of  the  vortex  sheet  originating  at  the  point  of  flow  separation.  The 
strength  of  the  vortex  is  also  related  to  the  orientation  of  the  vortex  sheet 
and  to  the  level  of  vorticity  shed  from  the  boundary  layer.  These  viscous 
phenomena  are  critically  affected  by  the  design  of  the  leading  edge. 

Obviously,  the  best  leading  edge  design  for  supersonic  cruise  may  not 
be  the  best  for  efficient  subsonic  or  transonic  speeds.  It  may  not  be  the 
best  for  transient  and  stability  effects.  Therefore,  a rounded  leading  edge 
may  have  to  be  designed  to  give  the  best  all-around  performance  and 
stability.  With  the  tools  to  analyze  the  rounded  leading  edge  in  the 
presence  of  a leading  edge  vortex,  the  most  efficient  all-around  leading  edge 
may  be  designed.  It  may  well  be  that  neither  stability  nor  performance  need 
be  compromised. 

The  available  theoretical  tools  for  wing/vortex  analyses  are  very  sensi- 
tive to  the  starting  conditions  at  the  leading  edge.  In  the  thin  wing 
potential  flow  solutions,  such  as  those  of  References  5 and  6,  it  is  correctly 
assumed  that  the  vortical  elements  and  velocity  vectors  at  the  wing  leading 
edge  lie  in  the  plane  of  the  wing,  but  their  orientation  in  the  plane  must  be 
assumed.  Techniques  using  quadrilateral  panels  for  the  feeding  sheet,  such 
as  that  of  Weber,  et  al.\  apparently  avoid  this  orientation  problem  for  thin 
wings.  For  thick  wings,  however,  the  panel  would  have  to  be  oriented  with 
the  actual  direction  of  the  feeding  sheet  leaving  the  boundary  layer.  The 

"Viscous  Box"  method  has  also  shown  (Sections  2.2  and  2.4)  that  the  orienta- 

tion and  strength  of  the  vorticity  coming  from  the  boundary  layer  into  the 
box  is  critical.  Therefore,  it  is  quite  conclusive  that  a boundary  layer 
analysis  at  the  wing  leading  edge  is  necessary. 

Additional  computations  must  be  made  to  represent  the  boundary  layer 
region  forward  of  the  box.  Strictly  speaking,  it  is  required  to  compute  a 
boundary  "region",  i.e.  a region  including  higher-order  effects  than  those 
occurring  in  first-order  boundary  layer  theory  (see  Figure  33).  In  order  to 
establish  with  accuracy  the  forward  edge  velocity  and  vorticity  values  for 

the  box,  it  is  necessary  to  perform  a computation  from  the  wing  attachment 
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Figure  33  Section  Showing  First  and  Second  Order  Boundary  Layer  Domains 


line  through  the  boundary  layer  and  into  the  second  order,  rapidly  thickening, 
boundary  layer  region. 

Recognizing  the  sensitivity  of  the  viscous  box  vortex  flow  region  to 
the  vorticity  imposed  by  the  leading  edge  boundary  layer  region,  it  is 
essential  that  three  dimensionality  be  accounted  for  in  any  model  of  the 
boundary  region.  For  this  purpose,  it  is  deemed  necessary  and  sufficient  to 
apply,  the  "infinite  yawed  wing"  assumption.  This  method  has  proven  quite 
fruitful  in  the  analysis  of  boundary  layers  near  attachment  lines  and  even 
for  some  complete  wings  of  sufficiently  high  aspect  ratio^3. 


This  boundary  layer  problem  will  be  investigated  in  follow-on  efforts  of 
this  program,  and  it  will  be  applied  to  the  thick  leading  edge  configuration. 
Only  the  thin  sharp  edge  wing  is  considered  in  this  report. 


The  thin  wing  is  not  without  leading  edge  problems  for  it  is  also 
necessary  to  know  the  strength  of  the  vorticity  being  shed  from  the  leading 
edge.  It  is  also  necesssary  to  locate  it  appropriately  on  the  leading  edge 
of  the  vorticity  box  to  assure  that  the  vorticity  is  shed  into  the  box.  In 
the  discussion  of  "Results  of  First  Iteration  for  a Thin  Delta  Wing," 
Section  3-8,  some  problems  related  to  this  input  vorticity  are  discussed. 


Obtaining  the  leading  edge  vorticity  for  the  thin  sharp  edge  wing  is  not 
as  readily  available  from  theory  as  it  is  for  the  thick  leading  edge.  Some 
theoretical  or  empirical  techniques  must  be  used  other  than  vortex  lattice 
to  determine  the  pressure  distribution  on  the  lower  surface  of  the  thin  air- 
foil which  in  turn  is  used  to  compute  the  boundary  layer  vorticity.  At  this 
time  the  theoretical  method  is  not  available,  and  vorticity  deduced  from 
laser  velocimeter  data  (see  Section  A.O)  is  used  for  input  to  the  thin  delta 
wing  investigation. 


3 . 7 Distributed  Circulation  from  Viscous  Box  Results 

After  the  first  iteration  of  the  hybrid  viscous  vorticity  box  method  and 
potential  flow  vortex  lattice  method,  circulation  distribution  can  be  defined 
from  the  vorticity  box  results.  These  distributions  can  in  turn  be  used  for 
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the  second  and  following  iterations  of  the  vortex  lattice  method.  A distri- 
buted free  vortex  system  was  discussed  in  Section  3.** 

Figure  3**  shows  the  lines  of  constant  vorticity  for  one  of  the  cases 
run  for  box  one  of  Figure  19-  The  isovars  shown  are  for  the  downstream  end 
of  box  one  or  the  upstream  face  of  box  two.  Positive  and  negative  discrete 
circulations  can  be  formed  as  shown  by  the  rotation  symbols  of  Figure  3^, 
where'  the  circulation  is 


r 


£ dA 


(18) 


and  A is  area  chosen  to  represent  the  vortex  core  circulation  or  portions  of 
the  vortex  sheet.  The  centroid  of  each  circulation  is  obtained  from  moments 
In  the  case  shown,  four  elements  make  up  the  vortex  sheet  and  one  element  for 
each  of  the  remaining  positive  and  negative  vorticities  make  up  all  other 
circulation.  The  vorticity  at  the  lower  right  corner  of  the  box  is  not  too 
realistic  because  this  part  of  the  box  is  too  close  to  a bound  vortex  element 
of  the  vortex  lattice.  This  gives  unrealistic  vertical  velocities  as  shown 
in  Figure  31*- 

It  is  noted  that  the  vortex  sheet  is  of  unusual  shape.  This  is  not  as 
expected  and  recommendations  for  correcting  this  are  given  in  Section  3-8. 

In  order  to  maintain  a realistic  representation  of  the  distributed  free 
vortex  system,  the  elements  chosen  to  make  up  the  vortex  sheet  of  box  one 
will  have  to  remain  constant  as  they  move  through  the  boxes.  It  will  there- 
fore be  necessary  to  investigate  the  vortex  sheet  at  each  box  face  to  obtain 
its  total  circulation,  then  divide  it  into  areas  that  will  integrate  to  give 
the  same  circulation  elements  as  previously  generated  in  the  preceding  face; 
so,  each  new  box  face  will  contain  each  of  the  elements  previously  generated 
plus  one  more  for  new  vorticity  being  fed  into  the  box.  This  method  will  be 
compatible  with  the  distributed  free  vortex  lattice  discussed  in  Section 
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3.8  Results  of  First  Iteration  for  Thin  Delta  Wing 

The  various  programs  and  the  modeling  techniques  for  integrating  the 
programs  have  been  given  in  Section  3-1  through  3-7-  The  modeling  technique 
has  been  applied  to  the  thin  delta  wing  planform  of  Figure  19.  Only  one 
angle  of  attack  (22°)  is  considered.  Both  the  single  and  distributed  free- 
vortex  systems  discussed  in  Section  3.1»  have  been  used  to  determine  the 
boundary  conditions  for  each  of  the  six  viscous  vorticity  boxes.  Viscous 
flow  computations  have  been  completed  for  all  the  boxes  in  both  cases. 

3-8.1  Input  for  Free-Vortex  Strength  and  Position  and  Leading  Edge  Vorticity 

To  begin  the  iterations  the  free  vortex  strengths  and  positions  are 
estimated  by  the  method  presented  in  Section  3*3-  Figure  23  shows  the 
estimated  free-vortex  circulation  strength  distribution  as  a function  of  the 
fraction  of  the  root  chord  Cr.  The  value  of  a = a/tany  * For  this 

level  the  values  of  Y/S  and  1/S  are  0.66  and  0.22,  respectively.  These 
positions  are  not  fixed  percentages  on  the  actual  case,  but  reasonable 
estimates  are  made  for  initiaiton  of  the  program.  It  was  discovered  after 
much  of  the  effort  was  already  accomplished  that  an  error  had  been  made  in 
reading  the  value  of  Y/S  at  .75  instead  of  0.66.  This  makes  the  vortex  too 
close  to  the  leading  edge  and  both  the  laser  data  and  vorticity  box  results 
indicate  the  vortex  core  is  further  inboard.  This  indication  from  the 
vorticity  box  is  a good  sign  that  convergence  is  possible. 

Another  input  required  is  the  vorticity  at  the  leading  edge  of  each 
vorticity  box.  This  information  was  not  available  from  theory  but  the  laser 
data  was  at  least  a good  approximation.  Laser  data  is  available  for  this 
estimate  from  only  one  box,  i.e.  at  the  face  of  box  2.  At  this  station,  the 
vertical  velocity  was  measured  along  lines  normal  to  the  wing  leading  edge  at 
various  heights  above  and  very  near  the  wing  surface.  Due  to  the  physical 
arrangement  of  the  laser,  measurements  could  not  be  made  right  at  the  surface 
level.  The  results  are  plotted  in  Figure  35.  If  points  could  be  obtained 
across  the  boundary  layer  on  a line  in  the  plane  of  the  wing,  the  vorticity 
would  simply  be  the  boundary  layer  edge  velocity  divided  by  the  boundary 
layer  thickness.  The  data  is  very  limited  to  accomplish  this,  but  a line 
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has  been  drawn  to  represent  crossing  the  boundary  layer.  The  slope  of  this 
line  is  the  vorticity,  at  least  the  major  component  of  vorticity. 

C(rp-)  = 3(V/Uj/3(n/c)  s 70  ( 1 9) 

^oo 

where  c is  the  width  of  box  2, 

n is  the  boundary  layer  thickness, 

V is  the  vertical  velocity. 

Variation  of  the  leading  edge  vorticity  along  the  leading  edge  of  the 
wing  is  estimated  from  the  logic  that  the  vorticity  is  directly  proportional 
to  the  circulation  growth  at  the  leading  edge  as  determined  from  vortex 
lattice  method.  The  values  of  the  growth  in  circulation  of  the  leading  edge 
vortex  element  are  given  as  a function  of  the  box  positions  in  the  following 
table.  Proportioned  vorticity  values  of  £ are  presented  also. 


TABLE  1 ESTIMATED  LEADING  EDGE  VORTICITY 


X/Cr 

ar/ u„ 

£/U» 

.2 

.486 

50.9 

-3 

.438 

45.5 

.382 

.358 

37.1 

-5 

.252 

26.1 

.578 

.195 

20.3 

-75 

.160 

1 6. 6 

.878 

.180 

18.7 

3-8.2  Resulting  Vorticity  Box  Circulation,  Velocity,  and 
Vorticity  Distribution 

Figures  38  through  38  shows  the  vorticity  contours  and  velocity  vectors 
as  computed  by  the  vorticity  box  method  for  the  aft  faces  of  boxes  1 through 
3.  The  results  are  not  satisfactory  at  this  time,  but  they  do  indicate  what 
is  necessary  to  make  them  satisfactory.  Boxes  4,  5,  and  6 results  are  not 
instructive  and  are  not  presented.  Difficulties  and  changes  required  are 
presented  in  the  following  discussion. 
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Beginning  with  box  one,  the  vorticity  contours  are  presented  in  Figure 
36.  The  results  here  are  qualitatively  very  good.  A feeding  sheet  can  be 
noted  as  well  as  a primary  vortex.  A relatively  strong  field  of  negative 
vorticity  can  also  be  found.  Centroids  and  integrated  strengths  of  a primary 
and  secondary  vortex  circulation  can  be  determined.  Also,  the  vortex  sheet 
can  be  divided  into  discrete  circulation  of  vortex  elements  for  input  into 
the  second  vortex  lattice  interaction. 

Quantitatively,  the  results  are  not  good.  From  an  integration  of  the 
vorticity  to  get  circulation,  the  total  circulation  is  1.202  with 
the  centroid  at  Y/C  * 0.275  and  Z/C  = 0.551,  where  c is  the  width  of 
the  box  and  Y is  from  leading  edge.  The  input  values  of  the  free-vortex 
centroid  were  Y/C  = .34  and  Z/C  =.30.  This  is  the  centroid  of  the  primary 
circulation  only.  The  primary  vortex  in  the  vorticity  box  is  at  approxi- 
mately Y/C  = .42  and  Z/C  = .23.  The  secondary  vortex  centroid  is  at  about 
Y/C  = .2  and  Z/C*.08.  As  noted  in  the  previous  section,  the  original 
estimate  of  Y/C  should  have  been  further  from  the  leading  edge.  It  is 
encouraging  that  the  vorticity  box  moves  the  primary  vortex  inboard. 

The  amount  of  circulation  from  the  vorticity  box  is  low  compared  to  the 
initial  input  which  was  T/U00  = 4.85.  r/U^  of  1.202  from  the  vorti.city  box  is 
the  total  circulation.  If  the  negative  circulation  is  subtracted  from  the 
total,  the  positive  circulation  is  approximately  1.488.  This  is  still  very 
low  compared  to  the  input. 


There  are  at  least  two  reasons  for  the  above  anomaly.  First,  and  most 
important,  is  described  by  referring  to  the  following  sketch. 

1 


The  vorticity  box  used  has  its  vertical  edge  at  the  leading  edge;  however, 
the  vorticity  contours  of  the  computed  results  indicate  the  vortex  sheet  is 
forced  to  have  a reverse  curve,  not  as  the  sketch  shows  it.  The  velocity 
vectors  on  the  vertical  edge  of  the  box,  figure  36,  show  very  little  inflow 
into  the  box;  therefore,  very  little  of  the  input  vorticity,  at  the  first 
vertical  mesh  point,  is  fed  into  the  box.  The  remedy  for  this  is  to  extend 
the  box  ahead  of  the  leading  edge  as  shown  by  the  dotted  lines  in  the  sketch. 
Now,  .the  input  leading  edge  vorticity  will  be  put  on  the  bottom  of  the  box 
just  ahead  of  the  leading  edge.  The  strong  vertical  components  of  velocity 
there  will  feed  the  vorticity  into  the  box,  and  the  correct  shape  of  the 
vortex  sheet  should  be  made.  This  proposed  correction  to  the  box  boundaries 
is  not  available  at  this  time.  It  will  probably  not  be  necessary  for  the 
thick  rounded  leading  edges  since  the  boundary  layer  turns  over  the  airfoil 
before  separating  except  at  very  high  angles  of  attack. 

The  second  reason  for  losing  vorticity  may  be  due  to  the  relatively 
large  mesh  used  for  these  trial  investigations.  The  effect  of  numerical 
accuracy  cannot  be  determined  until  the  problem  above  is  explored. 

Vorticity  contours  and  velocity  vectors  for  boxes  two  and  three  are 
shown  in  Figures  37  and  38.  These  suffer  even  worse  than  box  one  because 
of  the  same  reasons  as  well  as  the  fact  that  each  are  dependent  on  the 
preceding  box  results.  In  order  to  check  our  theory  that  vorticity  is  not 
being  fed  into  the  boxes,  the  leading  edge  input  vorticity  for  box  two  was 
moved  to  a position  where  strong  inflow  is  occurring  on  the  leading  edge, 
i.e.  at  vertical  grid  positions  6 and  7-  This  indeed  shows  a large  increase 
in  circulation,  about  three  times.  Figure  39  shows  the  vorticity  and 
velocity  vectors  for  this  case. 

It  was  decided  that  another  check  was  necessary.  Vorticity  entering 
a box  at  its  face  must  travel  downstream  on  the  order  of  3 or  k box  widths 
before  entering  the  primary  vortex  core.  Results  of  Figure  39  are  for  a 
downstream  distance  of  less  than  one  box  width;  so  this  case  was  run  again 
keeping  all  boundary  conditions  the  same  while  stepping  further  down  the 
box.  The  results  do  indicate  a rolling  up  of  vorticity  into  primary  vortex. 


65 


o u. 
— o 


Figure  39  Downstream  Box  2 Contours  and  Velocities 

with  Leading  Edge  Vorticity  Shifted  Vertically 
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The  use  of  rectangular  boxes  leads  to  difficulty  in  sustaining  a smooth 
numerical  computation  of  the  flow.  In  fact,  the  present  demonstration  indi- 
cates a breakdown  of  the  vortex  flow  pattern  following  transition  to  a larger 
box,  and  the  requirement  of  solution  interpolation,  as  discussed  in  Section 
3.5,  appears  to  be  a major  cause  of  this  behavior.  As  noted  earlier,  curvi- 
linear coordinate  system  is  needed  for  the  viscous  computation  in  order  to 
take  account  of  vortex  growth.  A spherical  system  would  probably  be  required, 
in  w^ich  case  terms  involving  two  curvatures  would  be  added  to  the  governing 
equations.  Such  a curvilinear  mesh  could  also  be  extended  forward  of  the 
leading  edge,  as  depicted  in  Figure  kO.  Thus,  a system  possessing  two 
curvatures  would  increase  in  physical  dimension  with  downstream  distance 
while  retaining  a fixed  number  of  mesh  points,  and,  most  importantly,  no 
discontinuities  would  be  encountered  in  the  stepping  direction. 

3.8.3  Pressures  from  Vorticity  Box 

The  most  realistic  results  to  date  have  been  obtained  from  box  one  of 
the  thin  delta  wing  investigation.  It  is  believed  the  pressure  calculation 
technique  as  described  in  Section  2.5  is  satisfactory  (at  least  for  the  fine 
grid  solutions).  Results  from  box  one  are  presented  in  Figure  1*1.  Absolute 
values  are  not  correct,  but  the  shape  is  typical  of  delta  wing  pressure 
distribution  except  .at  the  leading  edge.  Better  definition  of  leading  edge 
input  conditions  should  improve  the  leading  edge  pressures. 


A.  LASER  VELOCIMETER  EXPERIMENTS 


A. 1 General  Description 

Laser  velocimeter,  LV,  experiments  were  conducted,  primarily  to  determine 
the  velocity  vectors  for  the  boundary  conditions  at  the  face  of  one  of  the 
boxes.  The  delta  wing  and  LV  optics  are  shown  in  Figure  A2.  In  this  photo- 
graph the  LV  is  in  the  forward-scatter  mode.  With  the  wing  mounted  on  the 
side  wall,  the  laser  beam  comes  through  the  wind  tunnel  bottom  window  and  runs 
parallel  to  the  wing  upper  surface.  The  receiving  optics  are  looking  through 
the  side  window  from  a A5  angle  to  the  vertical.  In  this  mode  velocities 
parallel  to  the  tunnel  axis  and  normal  to  the  side  walls  are  obtained.  The 
LV  was  later  changed  to  back-scatter  mode  where  the  beam  comes  through  the 
side  window  and  the  receiving  optics  are  looking  through  the  same  window,  but 
off  slightly  to  the  side  of  the  laser  beam.  In  this  mode,  the  third  component 
of  velocity  (normal  to  tunnel  bottom  window)  and  again  the  component  parallel 
to  the  tunnel  axis  are  measured. 

A. 2 The  Model,  Wind  Tunnel,  and  Seeding 

The  model  is  a simple  flat  plate  delta  wing  with  65°  leading  edge  sweep. 
Leading  edges  are  sharpened  by  beveling  the  1/8"  (.3175  cm)  plate  on  the  lower 
surface.  A strut  from  the  model  lower  surface  to  the  side  wall  supports  the 
model.  Additional  support  is  provided  with  wires  off  the  lower  surface  to  the 
top  and  bottom  of  the  tunnel.  Root  chord  of  the  model  is  8.57  inches 
(21.768  cm). 

The  wind  tunnel  is  an  induction  tunnel  with  a centrifugal  blower  driving 
it  at  its  downstream  end.  Maximum  test  section  speed  is  approximately  120 
ft/sec.  (36.58  m/sec).  All  tests  have  been  run  at  50  ft/sec.  (15-2A  m/sec). 
The  test  section  is  approximately  29x17  inches  (73-66  x A3. 18  cm)  with  the 
large  dimension  in  the  vertical  direction.  This  facility  was  designed  and 
fabricated  primarily  to  develop  operating  techniques  for  the  LV  system. 

One  of  the  systems  under  study  is  the  method  of  seeding  the  tunnel.  In 
this  case  a Sommenist  ultra-sonic  atomizing  nozzle  is  used  upstream  in  the 
contraction  section  of  the  tunnel.  Water  seemed  to  give  good  seeding  rates 
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and  it  was  used  as  the  seeding  material.  Particle  size  has  not  yet  been 
determined.  This  will  be  done  once  the  purchased  particle  analyzer  arrives. 


k. 3 Results  of  Delta  King  Test 

Boundary  conditions  were  desired  at  one  box  to  obtain  typical  input 
conditions  for  the  viscous  vorticity  box.  Specifically,  it  was  desired  to 
measure  the  edge  velocities  on  the  box  face  and  velocities  on  a line  normal 
to  the  leading  edge  in  the  plane  of  the  upper  surface  of  the  delta  wing  and 
extending  from  the  leading  edge  out  to  free-stream  conditions.  Since  the  box 
face  is  normal  to  the  leading  edge,  the  single  line  ahead  of  the  wing  is  also 
in  the  plane  of  the  box  face. 

The  face  of  box  2 is  chosen  for  the  LV  analysis.  A sketch  of  the  model 
and  results  of  the  velocities  measured  along  the  extended  line  are  given  in 
Figure  35*  It  was  very  difficult  to  measure  exactly  on  a line  in  the  plane 
of  the  wing  upper  surface;  therefore,  the  velocities  are  given  for  several 
heights.  As  a first  approximation,  results  from  the  vertical  velocity 
distribution  across  the  boundary  layer  at  6/c=.07  provides  an  estimate  of  the 
leading  edge  vorticity;  i.e. 


5(c/Uj  = »(v/Uj/8(Y/e)  = 70 


(20) 


where  c is  the  width  of  the  box  face  (leading  edge  to  plane  of  symmetry). 
(Note,  for  investigations  presented  in  Section  2.0,  the  reference  chord  is 
.5  c;  therefore,  5/11^  = 35).  Values  of  U and  V are  based  on  tunnel  wind  axis 
coordinates,  and  it  would  be  slightly  more  correct  to  convert  to  model  box 
coordinates;  i.e.  U,  V,  W which  are  calculated  from  U',  V,  W by  equation 
(15)  of  Section  3-5.  The  data  of  Figure  35  have  not  been  converted  partly 
because  the  W component  is  not  available  and  partly  because  the  U component 
in  the  edge  boundary  layer  could  not  be  measured  at  6/c  = 0. 

The  next  objective  was  to  measure  the  boundary  conditions  around  the 
face  of  box  2.  This  has  been  done  twice.  The  first  time  only  the  U1  and  V' 
velocity  components  in  the  tunnel  coordinate  system  were  measured.  The  W* 
component  was  measured  at  a later  date  with  the  LV  in  the  backscatter  mode. 
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In  the  first  test  the  model  support  and  the  tunnel  seeding  methods  were  dif- 
ferent than  when  the  backscatter  was  done  on  the  second  test.  It  was  felt 
the  support  and  seeding  was  much  better  in  the  second  test;  so  forward 
scatter  (U1  and  V1  components)  were  measured  again. 

During  the  forward-scatter  measurements  of  the  second  test,  it  proved  to 
be  more  accurate  and  faster  to  take  data  in  planes  normal  to  the  tunnel  axis. 
This  is  faster  data  acquisition,  but  causes  much  more  analyses.  Three  planes 
of  data  grouped  in  the  area  of  the  box  coordinate  plane  were  taken  and  then 
cross  plots  are  made  to  get  all  data  in  the  box  coordinate  system  plane.  (The 
box  coordinate  system  plane  is  normal  to  the  wing  surface  and  wing  leading 
edge.)  Both  forward  and  back-scatter  data  were  taken,  i.e.  U1,  V1,  and  W 
components,  and  enough  data  has  been  taken  to  map  the  fiow-field  vectors  over 
most  of  the  box  face.  Only  the  edge  velocities  are  obtained  at  this  time  due 
to  time-consuming  data  analyses.  These  are  shown  in  Figure  43  with  the 
vectors  rotated  to  box  coordinates.  Figures  27  and  28  show  the  measured 
vectors  in  the  wind  axis  system  comparing  with  the  results  from  the  vortex 
lattice.  The  latter  show  the  initial  assumed  vortex  location  is  too  close  to 
the  surface  and  the  wing  leading  edge.  The  experimental  velocities  along  the 
leading  edge  of  the  box  also  indicate  the  vortex  feeding  shed  crosses  at 
about  0.25  Y/c  as  is  also  shown  in  Figure  43. 

Figure  43  also  shows  the  effective  vortex  center  as  obtained  from  LV 
data  by  finding  the  intersection  of  the  lines  for  zero  spanwise  and  vertical 
velocity.  A further  observation  of  the  LV  data  indicates  the  vortex  core 
extends  beyond  the  upper  bounds  of  the  assumed  box.  This  of  course  violates 
the  theory  and  will  have  to  be  corrected  on  subsequent  efforts. 
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5.  CONCLUSIONS 


1.  The  computational  viscous  flow  model  has  been  developed  to  provide  a 
qual i tat ively  valid  solution  to  the  vortex/wing  interaction  problem. 

2.  It  is  shown  how  the  viscous  method  can  be  combined  with  a potential 
flow  model  in  an  iterative  technique  to  produce  the  flow  field  for 
a delta  wing  at  high  angle  of  attack.  Some  difficulties  in  getting 
the  appropriate  boundary  and  starting  conditions  are  apparent,  but 
none  discredit  the  feasibility  of  the  technique. 

3.  Accuracy  of  the  viscous  box  method  appears  satisfactory,  but  as  of 
now  there  is  no  exact  solution  to  compare  with.  The  only  proof 
will  appear  when  results  of  a converged  solution  for  a finite  wing 
are  compared  with  accurate  experimental  data. 
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6.  RECOMMENDATIONS  FOR  FURTHER  EFFORT 

There  are  some  additional  efforts  already  planned  towards  the  fruition 
of  this  program.  These  are:  (l)  Provide  a leading  edge  boundary  layer 
program  to  compute  the  input  vorticity  and  boundary  layer  velocities  in  order 
to  make  the  program  applicable  to  thick  leading  edge  configurations,  (2) 
interface  the  boundary  layer  and  viscous  box  theory  with  a thick  wing 
potential  flow  theory,  and  (3)  include  the  secondary  vortex  in  the  potential 
flow  theory. 

There  are  further  efforts  not  yet  in  the  planning  stage  that  should  be 
accomplished  soon.  Among  these  are:  (l)  change  coordinate  system  to 
spherical  to  eliminate  stepping  the  box  size,  (2)  modify  program  to  allow 
the  viscous  box  to  extend  beyond  the  leading  edge  to  improve  the  thin  wing 
calculations  and  to  eliminte  high  angle  of  attack  limitations  on  the  thick 
wing  problem. 

Future  efforts  should  include:  (l)  the  complete  integration  into  one 
program,  all  of  the  related  viscous  and  potential  flow  models  described  above, 
(2)  the  method  should  be  applied  to  wings  other  than  deltas  and  to  other 
related  viscous  problems  such  as  "Upper  Surface  Blowing,"  Spanwise  Blowing, 
and  Afterbody  Drag,  (3)  a turbulence  model  should  be  integrated  into  the 
viscous  model,  and  (i*)  the  viscous  box  model  should  also  be  extended  to 
include  two  sides  with  no-slip  flow. 
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LIST  OF  SYMBOLS 


AR 

F 

H 

u,v,w 
u,v,w 
u1 ,v' ,w' 
c 


Cr 

S 


x,y,z 

5,n,C 


a 

r 

Y 

5 

v 

P 

u> 

(1 

A 


aspect  ratio 

column  vector  of  vorticity  components 
total  pressure  function 

velocities  in  the  x,y,z  directions,  respectively 
mean  velocities  in  the  x,y,z  directions,  respectively 
wind  tunnel  velocities  in  the  x,y,z  directions,  respectively 
box  chord,  normal  to  leading  edge 
wing  root  chord 

wing  span,  normal  to  plane  of  symmetry 
lift  coefficient,  vortex  lift 
pressure  coefficient 

streamwise,  vertical,  and  spanwise  coordinates 
vorticity  in  x,y,z  directions,  respectively 
angle  of  attack  of  wing 
ci rculation 

included  angle  between  wing  leading  edge  and  C|_ 

distance  normal  to  wing  surface  and  leading  edge 
boundary  layer  thickness 

kinematic  viscosity 

density  of  flow  in  wind  tunnel 

vorticity  vector 

denotes  vector  unless  otherwise  noted 
wing  leading  edge  sweep 


AND  TYPICAL  OUTPUT 


— _w.  lfc.„  i n 
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1 

1 

A.  General  Comments  on  the  Capabilities  of  the  Program 

The  program  solves  approximated  steady-state  Navier-Stokes  equations 
using  finite  differences.  Vorticity  and  velocity  vectors  are  employed  as 
dependent  variables.  The  vorticity  transport  equation  is  parabolized  in 
one  of  the  coordinate  directions  by  assuming  diffusion  to  be  neglibible  in  i 

this  direction.  This  direction  is  labelled  as  X in  the  programs  coding. 

With  .initial  values  of  vorticity  and  velocity  specified  at  x“0  plane,  the 
solutions  are  obtained  by  marching  forward  in  the  X-direction.  The  domain 
of  integration  is  box-shaped  and  one  face  of  the  box  (y  * 0 plane  in  program 
coding)  has  no-slip  boundary  conditions  applied  on  it.  On  the  other  faces 
of  the  box  vorticity  is  specified  as  zero.  The  box  size  must  be  sufficiently 
large  to  meet  this  condition.  The  velocity  values  on  the  edges  of  the  box 
are  computed  using  Green's  function  technique.  Alternatively,  these  values 

j 

can  be  specified,  for  example,  using  potential  flow  solutions.  The  program 
is  capable  of  numerically  simulating  vortex-type  flows  near  no-slip  wall  and 
flows  with  wall  jet  interaction. 

J ] 

Discussions  on  the  formulation  of  governing  equations  and  the  finite 
difference  equations  used  in  the  computations  will  be  found  in  References  1 
and  2.  In  the  following  sections  details  that  a user  should  know  for  running 
case  studies  with  this  program  are  presented. 

B.  References 


1.  R.  M.  Scruggs  and  C.  J.  Dixon,  "Vortex/Jet  Wing  Viscous  Interaction 
Theory  and  Analysis,"  Report  No.  LG76ER0007,  Lockheed -Georg ia  Company, 
Feb.,  1976. 

2.  C.  J.  Dixon  and  R.  M.  Scruggs,  "Further  Development  of  a Viscous 
Vortex/Wing  Interaction  Program,"  Interim  Report  to  Office  of  Naval 
Research,  Contract  N0001 4-7^“C-0151 , NR21 5~233 - 
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Compu t a t ( ona I Sequence 


INITIALIZE 

(Subroutines:  BEGIN,  RESTART) 


WRITE  VALUES  REQUIRED  FOR  CONTINUING 
COMPUTATIONS  ON  TAPE 


IF  STATION  X>0  COMPUTE  PRESSURE 
COEFFICIENTS  ON  SURFACE 


OUTPUT:  VORTICITY  AND  VELOCITY  COMPONENTS, 
PRESSURE  COEFFICIENTS  ON  SURFACE 


IS  DESIRED  X-STATION  REACHED? 


ADVANCE  X VALUE.  PREPARE  FOR  COMPUTATION 
AT  NEXT  X-STATION 
(Subroutines:  STORE,  RES1 , EXTNL) 


COMPUTE  VORTICITY  AT  INTERIOR  POINTS 
USING  ADI  METHOD 
(Subroutine:  VTCTY) 


COMPUTE  VORTICITY  AT  INTERIOR  POINTS 
USING  ADI  METHOD 
(Subroutine:  VELOC) 


COMPUTE  SURFACE  VORTICITY 
(Subroutine:  CORRC) 


Iterate  for  Convergence 


COMPUTATIONS  COMPLETED  FOR  CURRENT  STATION 
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D.  Input  Detai Is 


l 


Card 

Nr.  Variable  Significance  and  Format 

1 XC  Cumulative  value  of  x from  initial  or  starting  plane. 

Input  as  zero  to  start  computation.  Input  value  of  x 
up  to  which  computations  were  done  on  previous  sum  to 
restart. 

FORMAT:  Free  (MAIN  Program) 


2 


UINF 

VINF 

WINF 

REN 

Cl 

A 


(U 

(2) 

(3) 

(4) 

(5) 

(6) 


^00 

v» 

w* 

Reynolds  number 

av,  core  diameter  of  vortex  (p.  46,  Ref.  1) 
aj,  defines  jet  shape  (p.  46,  Ref.  1) 

Set  to  zero  if  no  jet  is  present. 


3 


Note:  If  jet  is  not  present, values  given  to  YZC1  and 
YZCZ  are  immaterial . A variable  EPA  when  set 
to  zero  ’tums-off'  the  image  vortex.  (Line  68, 
Subroutine  BEGIN. ) 

FORMAT  for  Card  3:  8F8.0  (Subroutine  BEGIN). 


YZA1 

y-coord i nate 

of 

YZA2 

z-coordinate 

of 

YZB1 

y-coord i nate 

of 

YZB2 

z-coordinate 

of 

YZC1 

y-coord i nate 

of 

YZC2 

z-coordinate 

of 

image  vortex 
image  vortex 
vortex 
vortex 

jet  axis  in  initial  plane 
jet  axis  in  initial  plane. 


4 


MM 

NN 

ITN 


IXN 


DX 

DY 

DZ 


Number  of  grid  points  in  y-di recti  on. 

Number  of  grid  points  in  z-direction. 

Maximum  number  of  iterations  for  vorticity 
and  velocity  (Recommended  value:  10  to  20) 

Number  of  x-stations  (including  the  starting  or 
restarting  station.  For  example,  to  advance 
computation  from  x = xj  to  x = Xj+Ax,  IXN  should 
be  2). 

Ax.  Since  the  procedure  is  implicit,  there  is 
no  restriction  on  step  size. 

Ay 

Az 

FORMAT  for  Card  4:  413,  E8.2,  2F8.0  (Sub.  BEGIN). 


Variable 


Significance  and  Format 


DYP 

DZP 

VORTLE 


Ay 1 . New  va 1 ue  of  Ay  when  box  s i ze  changes . 

(See  Fig.  19,  Ref.  2). 

Az1.  New  value  of  z when  box  size  changes. 
FORMAT  for  Card  5:  10F8.0  (Subroutine  RESTART) 

Leading  edge  value  of  vorticity  component  £. 
FORMAT:  F8.0  (Subroutine  RESTART). 

Potential  flow  velocity  values,  punched  in  the 
following  sequence. 


Component  up 

UP2,NN  uP3,NN  UPMM,NN 

UP2,NN-1  UP3,NN-1  UPMM,NN-1 

UP2,1  UP3,1 UPMM,1 

Component  vp 


(Sequence  as  above) 
Component  wp 


(Sequence  as  above) 


FORMAT:  10F8.0  (lines  18-*- 22,  Subroutine  RESTART). 
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E.  Definition  of  key  FORTRAN  VARIABLES 


ALPH 

SLAM 


ICON 


I 


( 

IMAP  ■{ 


CP(M,N) 

f(m,n,i) 

fo(m,n,i) 

fp(m,n,i) 

fa(m,n,i) 

V,  VO,  VP 
VA(M,N, I ) 

yz(m.k) 


Angle  of  attack,  a in  degrees  (defined  in  line  1A  of  subroutine 
RESTART) . 

Sweep  angle,  X in  degrees  (defined  In  line  15  of  subroutine 
RESTART) . 

= 1 Continuity  equation  will  be  used  to  compute  and 
over-write  u-component  of  velocity. 

= 0 Continuity  equation  will  not  be  used. 

(Lines  11,  12  MAIN.  Recommended  value  zero.) 

= 1 Green's  function  will  be  used  to  compute  boundary  velocity. 

= 0 Green's  function  will  not  be  used. 

Note:  If  Green’s  function  is  to  be  used , the  calling 
sequences  to  subroutines  must  additionally  be 
changed,  to  suppress  input  of  potential  flow 
velocity  values.  See  MAIN  program  in  listing. 

Pressure  coefficient  at  m,n. 

Components  of  vorticity  vector  at  current  x-station. 

Vorticity  components  at  previous  x-station. 

Vorticity  components  at  current  x-station  corresponding 
to  previous  iteration. 

Alternate  solution  of  vorticity  obtained  using  w = 0 on  the  wall. 
Used  in  determining  surface  vorticity. 

Velocity  - with  meanings  similar  to  F,  FO,  and  FP. 

Velocity  values  corresponding  to  FA. 

y and  z values  of  modal  points. 


F.  Tape  or  Mass  Storage  File  Requirement 

One  tape  or  mass  storage  file  with  a designated  name  TAPE10  is  required. 
Value  of  XC,  DY,  DZ,  V and  F are  written  on  this  file.  These  values  are 
used  for  restarting  the  computations  and  for  usage  subsequent  machine 
plotting. 


| 


86 


G.  Table  of  Subroutines  and  Their  Purpose 


Subroutine 

Name 

Purpose 

Other 

Subroutines 

Used 

1.  BEGIN 

Read  in  key  parameters  - Initialize 
vorticity  and  velocity  fields. 

None 

ENTRY-EXTNL 

Advance  x-station  and  re-ini tial ize. 

None 

2.  CONTN 

Compute  surface  vorticity. 

None 

3.  CORRC 

Compute  surface  vorticity. 

None 

k.  EQTNS 

Form  coefficients  in  finite  difference 
equations.  That  is,  matrices  [A], 

[B],  [C],  [D]  of  Eq.  (8),  Ref.  1. 

None 

5.  GUESS 

Make  an  initial  guess  for  velocity 
and  vorticity. 

(Not  used) . 

None 

6.  INVRT 

Invest  a 3 x 3 matrix. 

Augmented:  A Matrix  to  be  inverted 

3 Inverse 

IE  Test  parameter  which 
returns  a value  of  zero 
if  determinent  of  matrix 
is  less  than  10-15» 

None 

7.  OUTPT 

Points  out  computed  results. 

None 

8.  SOLVE 

Inverted  tri-diagonal  block  matrix 
using  Cholesky's  method. 

INVRT,  MULTY 

9.  MULTY 

Multiplies  two  matrices  (3X3)- 

None 

10.  STORE 

Sets  V0=V  and  F0  = F before  advancing 
x-station. 

None 

ENTRY-ST0RE2 

Sets  VP  = V and  FP  = F before  advancing 
i terat ion. 

11.  VELOC 


Computer  velocity  values  at  interior 
points  using  ADI  method. 


SOLVE 


12.  VTCTY 


Computes  vorticity  values  at  interior 
points  using  ADI  method. 


EQTNS 

MULTY 

SOLVE 
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Subroutine 

Name 

13-  DVEL 

14.  CPDIS 

15.  RESTART 

ENTRY-REST 

ENTRY-SHIFT 


Purpose 


Computes  velocity  on  boundaries  using 
Green's  function. 

Computes  pressure  distribution 
(Page  11-12,  Ref.  2). 

Reads  potential  flow  velocities  and 
converts  them  into  box  coordinates, 
interpolates  field  values  when  box 
size  changes. 

Reads  potential  flow  velocities  when 
there  are  intermediate  stations  within 
same  box. 

Writes  on  TAPE10  the  values  of  XC,  DY, 
DZ,  F and  V. 


Other 

Subroutines 

Used 


None 


None 


T 
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TkX  ■ .006751  crx  ■ 0.000000 

T«Z  > -.000067  crz  ■ 0.000000 
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